#### This script will plot microCLIM vs BIOCLIM
#### Establish directories

mc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/microCLIM/'
bc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/BC6'
out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/Plots/SurfaceComparison/'

#### List of surfaces

mc.files = list.files(mc.dir, recursive=T, full.names=T, pattern='.asc.gz')
bc.files = list.files(bc.dir, recursive=T, full.names=T, pattern='BC')[1:4]

# Names for surfaces

mc.names = gsub('_','',gsub('.asc.gz','',list.files(mc.dir, recursive=T, pattern='.asc.gz')))
bc.names = gsub('_','',gsub('.asc.gz','',list.files(bc.dir, recursive=T, pattern='BC')[1:4]))

### Cycle through surfaces

for (num in c(1,4,5,6))

{

### Plot daily Obs vs. Preds on same graph for max and min temp

t.mc = read.asc.gz(mc.files[which(substr(mc.names,nchar(mc.names),nchar(mc.names))==num)])
t.bc = read.asc.gz(bc.files[which(substr(bc.names,nchar(bc.names),nchar(bc.names))==num)])

t.mc = t.mc[which(is.na(t.mc)==F)]
t.bc = t.bc[which(is.na(t.bc)==F)]

png(paste(out.dir,'Surface 0',num,'_Comparison.png',sep=''), units='in', width = 5, height=5, res=500)

lm1 = lm(t.bc ~ t.mc)

newdata1 = data.frame(obs=t.mc) # Bind x-data for plots into a dataframe
	
rawpoly = predict(lm1, newdata1, interval="confidence"); rawpoly = cbind(newdata1,rawpoly) # Calculate 95% CI's for all x values in the AWAP linear model
	
polybottom = rawpoly[c(order(newdata1)),c(1,3)]
names(polybottom) = c('x','y')
polytop = rawpoly[c(rev(order(newdata1))),c(1,4)]
names(polytop) = c('x','y')
polyout = rbind(polybottom,polytop) # Calculate a polygon that represents the 95% CI's of the min linear model
		
lims = range(c(polyout[,1],polyout[,2])) # Calculate limits for x and y axes

plot(range(polyout[,1],na.rm=T),range(polyout[,2],na.rm=T),xlim = lims, ylim = lims, xlab = 'microCLIM (Observed)', ylab = 'BIOCLIM', main = paste('Surface 0',num,' Comparison',sep=''), type='n')

abline(b=1,a=0,col='black', lwd=2, lty=2) # 1:1 line
	
#polygon(polyout[,1],polyout[,2], col='#FF0000', lty=1, lwd=1, border=NA) # Plot a polygon representing the 95% CI's for the min linear model

#points(t.mc, t.bc, col='black')
		
abline(b=summary(lm1)[[4]][2],a=summary(lm1)[[4]][1],col='red',lwd=2, lty=1) # Plot a line representing the relationship between min and Empirical data		
	
dev.off()


}





